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Abstract 

We introduce a general characterization of sudden cosmological singu- 
larities and investigate the classical stability of homogeneous and isotropic 
J~ I . cosmological solutions of all curvatures containing these singularities to 

small scalar, vector, and tensor perturbations using gauge invariant per- 
turbation theory. We establish that sudden singularities at which the scale 
factor, expansion rate, and density are finite are stable except for a set 
^ ' of special parameter values. We also apply our analysis to the stability 

, of Big Rip singularities and find the conditions for their stability against 

(N . small scalar, vector, and tensor perturbations. 

PACS: 98.80.-k 
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I. INTRODUCTION 

o 

Q\ ' We have shown in [Ij that general relativistic Friedmann-Robertson- Walker 

■ (FRW) universes allow finite-time singularities to occur in which the scale factor, 

a(t), its time derivative, a, and the density, p, remain finite whilst a singularity 
occurs in the fluid pressure, p — ► +00, and the expansion acceleration, with 
a — + —00. Remarkably, the strong energy condition p + 3p > continues to 
hold. Analogous solutions are possible in which the singularity can occur only in 
arbitrarily high derivatives of a(t), [2J. This behaviour occurs independently of 
the 3-curvature of the universe and can prevent closed FRW universes that obey 
the strong energy condition from recollapsing [3] . These singularities can be seen 
in a wider context by classifying the behaviors of FRW universes containing 
matter with a pressure-density relation defined by p + p = "fp x , as shown in [J], 
and reviewed further in [5]. The sudden singular behaviour found for a range 
of values of (7, A) also encompasses the evolution found in a number of simple 
bulk viscous cosmologies studied in [6] . Subsequently, a number of studies have 
been carried out which generalise these results to different cosmologies, other 
density-pressure relationships, and theories of gravity [71(5]. 
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Other studies [5] have also related the sudden singularity behaviour to the 
formal classifications of 'weak' singularities according to the definitions intro- 
duced by Krolak [TU] and Tipler [TT], investigated the behaviour of geodesies, 
classified the other types of future singularity that can arise during the expan- 
sion of the Universe [T^], discussed the role of different energy conditions on 
different realizations of sudden singularities [121 [7] , and explored some obser- 
vational constraints on their possible future occurrence in our visible Universe 
|14j . The Weyl invariant will not diverge on approach to a sudden singularity 
(and there is no geodesic incompleteness [IS]), so it may represent part of a 
'soft' future boundary of the universe with low gravitational entropy - which 
could be as close as 8.7 Myr in the future. 

Most recently, the effects of quantum particle production have been studied 
and have been found to leave sudden singularities in place [16]. Specifically, 
it was shown that quantum particle production does not dominate over the 
classical background density on approach to a sudden singularity and does not 
stop it occurring or modify its properties, as can be the case for the Big Rip 
future singularities [T7J, [TH]. The effects of loop quantum gravity have been 
studied in cosmologies exhibiting classical sudden singularities and they may 
remove the sudden singularity under certain particular conditions [19] . Sudden 
singularities have also been studied due to their occurrence in various theories 
of modified gravity |20j . and we would only expect these modifications to be 
significant in this respect if they also dominate over general relativity effects at 
late times. 

In this paper, we will extend previous studies by investigating the classical 
stability of sudden singularities with respect to small inhomogeneous scalar, 
vector, and tensor perturbations using the gauge invariant formalism introduced 
by Mukhanov [3T] . We introduce a new characterization of sudden singularities 
in terms of the series expansion of the expansion scale factor on approach to 
the singularity. We show that, except for a subset of special parameter choices, 
they are stable to small perturbations if and only if the density does not diverge 
near the singularity. The latter is characteristic of sudden singularities. We also 
extend this analysis and apply it to 'Big Rip' singularities and determine the 
conditions under which they are stable and unstable. 

II. BACKGROUND THEORY 

Following [22], we consider perturbations of a general FRW metric, 

where K — 0, 1, —1 depending on whether the three-dimensional hypersurfaces 
of constant r time are spatially flat, closed or open. Here, r denotes conformal 
time, and is related to the comoving proper time, t, by a dr = dt. We shall also 
assume that the primary component of matter is a perfect fluid with energy- 
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momentum tensor (index notations run 1 < i, j < 3, < a, [3 < 3) 

T$ = (p+ P )u°>u -p6%, (1) 

but we shall not necessarily be specifying a specific equation of state linking p 
and p. 

For this background universe, the equations for the scale factor a(i), density 
p and pressure p are as follows (in units with 8ttG = 1 and c = 1), where the 
overdot denotes d/dt: 

p = -3-(p + p) (2) 



a 



^ = §"5 (3) 



(4) 



d ( p + 3p 



a V 6 



and we shall assume an equation of state with a functional form p = p(p). We 
can see by inspection that these equations permit finite-time singularities such 
that a, a and p are finite but p and d diverge as t — > t s . 

We shall adopt the following definition [53] : a sudden singularity will be said 
to occur at time t = t s if the scale factor a(t) can be written in the form 

a(t) = c Q + ci(t s - t) Xl + c 2 (t s - t) X2 +... (5) 

in a generalised power series about t s , where the Cj, are real constants, with 
Ci 7^ 0, co > and < Ai < A2 < . . . with at least one of the Ai non-integral 
(so that some derivative of a{t) blows up near the singularity). Note that the 
series ([5]) need not be infinite. Also, our form of §5§ assumes that the sudden 
singularity occurs in the future; however, the stability results we derive will also 
hold for past sudden singularities, defined in the analogous way using 

a(t) = c Q + ci(t - t s ) Xl + c 2 (t - t s ) X2 +... (6) 

This characterization encompasses the particular finite series expressions for 
FRW models with sudden singularities introduced in Refs. [UH], those arising 
in the solutions found in Ref. [4], and those studied in Refs. 

The energy-momentum tensor defined in (fTJ) leads to the following gauge- 
invariant perturbations [21] 

<5T ° = Sp, 5T° = -(Po+ Po)Su h ST; = -5p5) = -p\p)6p5). (7) 

Here, Sp, Sp and Sui are the gauge-invariant perturbations of the density, 
pressure and velocity. We shall now consider separately the behaviour of the 
tensor, vector and scalar modes as t — > t s . 

III. TENSOR PERTURBATIONS 
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Under tensor perturbations, the most general form of the line element is 

K , 



ds z 



cl\t) dr 



dx l dx 3 



where h\ = 0, h]y = 0, where a slash indicates a covariant derivative with respect 
to the spatial 3-metric. Note that the quantity hij is gauge-invariant. 
In conformal time, the equation for tensor perturbations is |24j 



2d 



2Khij = 0, 



(8) 



where a prime denotes a derivative with respect to conformal time. The cosmic 
time analogue of this equation is 



a ha + Zaahn + k hi 



2Kh lj = 



(9) 



for a plane wave perturbation with wavenumber k. We now set hij — veij, 
where is a time-independent polarisation tensor. This leads to the differential 
equation 

0. (10) 



a 2 v + 3adv + (k 2 + 2K)v 
From ([5]), we see that in the limit t — * t s , we have 

a{t) = c + Cl {t s - t) Xl + ... 



(11) 



a(t)=-c 1 \ 1 (t s -t) Xl - 1 + ... (12) 

where Ai > 0, and Cq,c\ ^ 0. We can substitute these into (|10p to obtain the 
asymptotic ordinary differential equation 

civ - ScoCiAiT^- 1 ?) + (fc 2 + 2K)v = 0, (13) 

where we have taken T = t s — t and dots now indicate differentiation with 
respect to T. We would like to investigate whether solutions to (fTB")) exhibit 
blow-up near T = 0. 

In the case Ai < 1, a very similar analysis to that described in the Appendix 
shows that v tends to a constant as T — > 0. Equations (|6T[) onwards still hold. 

In the case Ai > 1, all coefficients of (fl3|) are non-singular, so the ordinary 
differential equation is regular and hence has no singularity at T — > 0. In fact, 
if k 2 + 2K > 0, we end up with a simple harmonic oscillator (which is damped 
if Ai = 1), whose solutions are bounded at T = 0. 

Hence, both long- and short-wavelength tensor perturbations are bounded 
near the singularity. This clearly holds for all values of the A i; and the sudden 
singularity is always stable against inhomogeneous gravitational-wave pertur- 
bations. 



IV. VECTOR PERTURBATIONS 
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For vector perturbations, the most general form of the line element is 

K 



ds 2 = a 2 (r) dr 2 + 2S l dx l dr 



l + -( x ' + y * +z '))5 ij -F Mj 



dx l dx 3 



where S. 



1 1 . The quantity V% = Si — F[ is gauge- invariant. 



The only part of the energy-momentum tensor which contributes to vector 
perturbations is [5T] 



1 



& T i = -(Po+Po)Su ±l , 



(14) 



where Su±i is the part of Siii with zero divergence, and p and po are the 
background density and pressure, respectively. The equations for the vector 
perturbations are 



[2K - k 2 )V t = 2a(p +po)6u±_i 



(15) 



a(Vi,j + Vj. t y + 2a(V itj + V jti ) = 0. (16) 

For both long- and short-wavelength perturbations, we obtain as a conse- 
quence of the conservation of angular momentum: 

Vi — const x — , 

8v l = const x —t- — -, (18) 

where the physical velocities Sv z arc defined by Sv l = —a~ 1 5u±i. Hence, vector 
perturbations of the metric are bounded on approach to the sudden singularity 
since a — > a s < oo, p — > p s < oo and po — * oo. 

V. SCALAR PERTURBATIONS 
A. Overview 

We have shown, above, that tensor and vector perturbations do not diverge 
near the sudden singularity, for all values of in the form ([5]). The analysis of 
scalar perturbations is slightly more involved, and we shall need to consider var- 
ious cases according to the values taken by Ai and A2. This is not unexpected. 
The sudden singularity is primarily created by the behaviour of the pressure 
and so we expect the scalar perturbation modes associated with pressure inho- 
mogeneities to play a significant role in controlling the stability. 

Under scalar perturbations, the most general form of the line element is 

ds 2 = a 2 (r) ( (1 + 2(f>)dT 2 - 2B, l dx l dT - 
and we can define the usual gauge-invariant quantities 



(17) 



(1 - 2^0 ( 1 + ^-{x 2 + y 2 + z 2 ) ) 5ij + 2E. 



dx l dx 3 
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$ = 4> - - MB - E')}' and* = tp + —{B - E'). 
a a 

The equations for scalar perturbations are, following [22 : 

$ = * (19) 



a 2 $ + (4 + 3p'(p))aa$ + (2ad + (d 2 - K){\ + 3p'(p)) + p'(p)k 2 )<S> = (20) 

for plane wave perturbations with wavenumber fc, where p'(p) = dp/ dp. We can 
then find the gauge-invariant perturbed quantities as follows: 

6a- /2fc 2 + 6d 2 \ , x 

6p = p\p)6p. (22) 

B. The general case 

. First, consider the case when Ai = A is non-integral. In this case we find, in 
the limit t — > t s : 

a{t) = c + Cl (t s - t) x + ... (23) 

a{t) = -cx\(t s -tf' 1 + ... (24) 

o(t) = Cl \(\-l)(t s -t) x - 2 + ... (25) 
and from ([3]) and (HJ) we obtain the leading-order approximations 

P= 3 4 + 3 -^-(ts~t)^ + ... (26) 

c c 
Co C 

for the density and pressure, respectively. (Note that, regardless of the value of 
K, the density diverges if and only if A < 1, and the pressure diverges if and 
only if A < 2.) We also find that 

dp dp/dt 1 f d 3 — a 2 a + Kd\ cq(A — 2) x 



^ ^ dp dp/dt 3 \d(ad — d 2 — if) / 3ciA ^ s ^ ^ ^ 

to leading order, regardless of the value of K. 
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Substituting the above forms into (|2T))) . and neglecting higher-order terms, 
gives us the following equation (where T = t s — t and dots now indicate differ- 
entiation with respect to T): 

$ - (A - 2)T- 1 $ + ( — A 2 T A ~ 2 - ( fc2 - 3 ^)( A - 2 ) r -A _ $ = (2Q) 
V c o 3c ciA eg/ 

and we can neglect the last term since, for all values of A, it will be 

dominated in magnitude by one of the other two terms making up the coefficient 
of $, as T — » 0. We will now consider the cases A < 1 and A > 1 separately. 
If A < 1, then (2§J) becomes 

$-(A-2)T" 1 $ + — A 2 T A ~ 2 $ = (30) 
co 

and applying the substitutions P = T( 1 ~ A -'/ 2 <I> and x — CT X I 2 gives the equa- 
tion 

where primes denote differentiation with respect to x. If c\ > 0, we set C — 
2\/ c i/ c 0j and obtain a Bessel equation with solution 



P(x) = AJ u (x) + BY u (x) 
A- 1 

with A, i? arbitrary constants. Since we are considering the limit x — > 0, the 
leading-order solution is 

P(x) = Ax^r + Bx^r (32) 

or 

$(t) = A(t s -t) A_1 +B (33) 

where A, 2? are new constants. Therefore, since A < 1, $(f) diverges in general 
as t — > i s . The case ci < can be treated similarly. Note that this result holds, 
independently of k, for both long- and short- wavelength perturbations and for 
all K. 

Now suppose that A > 1. Then (f2"5)) becomes 

$ - (A - 2)T~ 1 $ - Z?T" A $ = (34) 
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where D = (fc2 2) . Analogously to the above, we substitute P = T (1 - X ^ 2 <i>, 

and then x = CT X ~ X / 2 , which leads to the Bessel equation 



(1-A) 2 
(2 -A) 

A-l 



P"x 2 + P'x+ [ x 2 -)^>)p = (35) 



2- A 



where we have set C = 2y / \D\/\X — 2|. 

As x 0, this has the asymptotic solution 

P(x) = Ax^i + Bi" (36) 

or 

= A(t s -i) A_1 +B. (37) 

Thus, if A > 1, the scalar metric perturbations do not diverge as t — > t s . 
Again, note that this holds for all values of K, and for both long- and short- 
wavelength perturbations. 

Hence, the scalar metric perturbations diverge for all < Ai < 1, and are 
bounded for all non-integral Ai > 1. In fact, the analysis above also holds for 
all integral Ai > 2. In this case, the density and pressure tend to constants near 
the singularity, and the metric perturbations do not diverge. Thus, we just need 
to deal separately with the boundary cases of Ai = 1 and Ai = 2. 

C. The case Ai = 1, A 2 ^ 2 

Consider the case Ai = 1 < A2 < . . . and assume first that A2 7^ 2. So, as 
t -» t s : 

a(t) = c + ci (t a - t) + c 2 (t s - t) A2 + . . . (38) 
a(t) = -c 1 -c 2 \ 2 (t s -t) X2 - 1 + ... (39) 
a(t) - c 2 A 2 (A 2 - l)(t s - t) x *- 2 + ... (40) 
a(t) = -c 2 A 2 (A 2 - 1)(A 2 - 2)(t s - t) x ^ 3 + ... (41) 

+ (42) 

c o 

p _ (K + 4) 2c 2 \ 2 (\ 2 -l)(t s -t) x *- 2 | ^ 

Co c 

We can see that p tends to a constant as we approach the singularity, but p 
tends to a constant if and only if A 2 > 2, and diverges otherwise. 



8 



Wc can also obtain the following expressions: 



P'(p) = 



ifK + cl^O, and 



3c 

CoC2A 2 (A2-l)(A2-2) u _ ,\A 3 -3 

&CgC2—cf—Kc 1 
3 Cl (K+cl) 

~ 3 



if 1 < A 2 < 2 

if 2 < A 2 < 3 

if A 2 = 3 

if A 2 > 3 



6c\ 



for all A 2 , if K + c\ = (i.e. K = -1, d = ±1). 

First, we consider the case K + c\ ^ 0. Substituting the forms of p'(p) into 
CP gives 

<i> + AT- 1 $ + BT- 1 <l> = if 1 < A 2 < 2 
< 5 + ^ T A 2 -3 ( j > + Br A 2 -3 $ = o if 2<A 2 <3 

$ + A$ + B$ = if. A 2 > 3 

where A and -B are always constants. The first differential equation yields an 
asymptotic solution <I>(t) A' + B'(t s — i) A2_1 , so the scalar perturbations 
are bounded in this case. Showing that solutions of the second differential 
equation are bounded as T — > is more cumbersome: a derivation is given 
in the Appendix. For A 2 > 3, the coefficients of the differential equation are 
constants, so the general solution has the form 

$(t) « A 1 exp(Ai(* s - t)) + B' exp(A 2 (i s - <)), (44) 

which is bounded as t — > t s . 

The same analysis can be used to show that the solution in the case K + cf = 
is also bounded as T — > 0. Notice that these results hold for both long- and 
short-wavelength perturbations. 

D. The case Ai = 1,A 2 = 2 

Now let Ai = 1, A 2 = 2 < A 3 < . . . As t -> t s : 

a(t) = c + ci(i s - t) + c 2 (t s - t) 2 + c 3 (t s - t) X3 + ... (45) 



a(t) = -ci - 2c 2 (t s - i) - c 3 A 3 (i s - f) 



A:j 



(46) 



o(*) = 2c 2 + c 3 A 3 (A3 - l)(t, - t) x "- 2 + ... 
a(t) = -c 3 A 3 (A 3 - 1)(A 3 - 2)(t s - i) A ^ 3 + . . 



P 



3(K + 4) 



(47) 
(48) 

(49) 
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, = (50) 
Co c 

so the density and pressure tend to constants near the singularity. Also: 

CoC3A 3 (A3-l)(A3-2) /■ f \A 3 -3 :f ^ \ s- 1 

3ci ( c ?+a:-2c 0C!1 ) " 2<A 3 <3 

J 1 3ci(c|+K— 2c c 2 ) J 

3(2 C0 c 1 2 -5-K) if A3 > 3 

and we can now substitute these into (f2"0")) to get the following functional forms 
of the scalar perturbation equation (where T = t s — t, dots indicate derivatives 
with respect to T, and A, B are constants): 

$ + AT X *- 3 <I> + BT X3 - 3 $ = if 2 < A 3 < 3 
$ + A$ + B$ = if A 3 > 3 

Equations of this form arose in the previous section, and we showed that their 
solutions are bounded as T — > 0. 

The above expressions for p'(p) are not well-defined if K = 2c$C2 — c\. For 
this special case, we can calculate 

'M_J - £S ^ a ( i »-*)" 1 if 2<A 3 <3 

and the corresponding scalar perturbation equation is 

$ + A$ + B$/T = 0. 

Solutions of this differential equation can be expressed in terms of Whittaker 
functions, and are bounded as T — > 0. 

Hence, the scalar perturbations are bounded for Ai = 1, A2 = 2. This is 
quite a strong result, since most 'nice' functions a(t) can be expressed as power 
series. We have shown that for all such functions, as long as the coefficient of the 
t term is non-zero near the singularity (i.e. a(t s ) ^ 0), the sudden singularity 
is a stable solution of the FRW equations. 

E. The case A x = 2 

The final case to consider is: 

o(t) = c + Cl (t, - tf + Oi% - t) x * + . . . (51) 

where the Cj 7^ 0, c > and 2 < A 2 < . . . with at least one of the A,; non- 
integral. As t — > t s , we have: 

d(i) - -2ci(i a - t) - c 2 X 2 (t s ~ t)^- 1 + ... (52) 
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o(t) = 2 Cl + c 2 A 2 (A 2 - l)(t s - t) x >- 2 + ... (53) 

o(t) = -c 2 A 2 (A 2 - 1)(A 2 - 2)(t, - i)^ 3 + . . . (54) 

and p -> 3K/c$,p — > -(4cic + K)/($. 

First, let us assume if = 0. There are three cases to consider for p'(p): 

/ _ coc 2 A 2 (A 2 -l)(A 2 -2) ^ _ t)A2 _4 j f 2 < A 2 < 6 

PV) = | 2c? 3" e g C2 fe-*) 2 if A 2 = 6 

I |^(*.-*) a if A 2 >6 

We now substitute these expressions into the equation for scalar perturbations 
(j2"0)) . Using A, .B, ... to denote constants, the functional forms of the differential 
equations obtained are: 

$ + AT X *- 3 <I> + (BT X *- 2 + Cfc 2 T A2 - 4 )$ = if 2 < A 2 < 4 
<£ + AT& + B$ = if A 2 > 4 

The solutions of the second differential equation are bounded as T — > 0. 
If we assume that k ^ 0, the solution of the first differential equation can be 
expressed in terms of hypergeometric functions, and is bounded as T — > 0. 

Now suppose that K ^ 0. There are two sub-cases: K ^ 2ciCo and If = 
2c\Cq. For the first sub-case, we have: 

C ^fK:^ 2) (t s -t) x ^ if 2<A 2 <4 
p'(p)={ 2i f K % 2Kc \ if A 2 = 4 

3(2c* -K) if A 2 > 4 

and the corresponding differential equations are 

q> + AT X2 - 3 <f> + (BT X2 - 2 + Cfc 2 T A2 - 4 )$ = if 2 < A 2 < 4 
$ + AT$ + B$ = if A 2 > 4 

with the same stability results as above. 

Finally, we treat the special case K = 2c\cq. We obtain: 

_BiI|^l(i s _f)-2 if 2 <A 2 <4 
-jz-{t s -t)- 2 if A 2 >4 




6c 1 C r _ 1 . + 3(V-^)C r _ 2$ = 

for all A 2 > 2, where p'(p) = CT~ 2 . In general, the solutions of this equation 
are of the form $(T) = T 7 where 
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7 2 + (-6 Cl C/c - 1) 7 + 3(fc 2 - K)C/c 2 = 0. 

For boundedness we need the two roots of this equation to be non-negative, 
i.e. we require —&c\C < cq and (k 2 — K)C > 0. It is easy to see that the first 
condition only holds for A 2 < 3. So, in this sub-case, we only have boundedness 
if A 2 < 3 and (fc 2 — 2c\Cq)C > 0, and divergence otherwise. 

So all the perturbations are bounded except for the very special case K 7^ 
0, K — 2ciCo, A2 > 3 (and perhaps A2 < 3). 

F. Consideration of the sign of c 2 

Note that our analysis here has not assumed a simple fluid equation of state 
of the form p = wp, or indeed any functional relation between the density and 
pressure of the matter source, as would characterize a perfect fluid of k-essence or 
its generalizations (see [25] for a discussion) . Thus there is no general constraint 
arising from the positivity of the square of an effective speed of sound, c 2 . In 
cases that reduce to fluids, or to k-essence and its relatives, it should be possible 
to introduce further constraints upon the series coefficients in order to preserve 
the positivity of c 2 . 

VI. EXTENSIONS TO BIG RIP MODELS 

We can extend the above analysis easily to Big Rip models [26], where the 
scale factor behaves as: 

a(t) = c (t s - i)"o + ci(t a - t)'" 1 + . . . (55) 

where rj < 0, % < r\ x < . . ., and also cq > 0, Cj ^ for all i. Then, to leading 
order, we have 

a(t) = -J] c (t s - t)"o-i + . . . (56) 

S(*) = VoiVo - l)co(t s - ty°- 2 + . . . (57) 

a(t) = -VoiVo - 1)(% - 2)c fe - t) r ' ~ 3 + • • ■ (58) 

P ='ii 1 l{t s -t)- 2 + ... (59) 

p= (2-3r 1o )r lo (t s -t)- 2 + ... (60) 

p'( P ) = 2/3% - 1 + . . . (61) 

Note that a(t), p,p —> 00 as t — > t s . It is easy to see that, near the singularity, 
the tensor perturbations hij and the vector perturbations Vi decay to zero. 
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However, the physical velocities 8v l are proportional to (t s — t)~ 2 ~ 4r * a , so for 
them to be finite as t — > t s we need r] < —1/2. 

Scalar perturbations obey the following equation: 



■■ . s $ f2(K-k 2 ) 2(k 2 -3K)\ 2 „ 



(62) 



where T = t s — t and dots denote differentiation with respect to T. A solution 
of (|B2")l can be found in terms of Bessel functions, and we find that the solutions 
asymptotically tend to 

$(T) = d +C 2 T- 1 - ,1 °. (63) 

So there is divergence as t — > t s if ry > —1, otherwise <£> is bounded. 

Note that these results match those of [37] , which shows that, for equations 
of state of the form p = ap (where a is a constant), there is a discrepancy in 
behaviour between the cases a > —5/3 and a < —5/3. It can easily be checked 
that these correspond to the cases rj < — 1 and ?7q > — 1 respectively. 

VII. CONCLUSION 



We have produced a simple general characterization of sudden singularities 
in FRW universes. By the use of gauge invariant perturbation theory we have 
investigated whether the existence of sudden singularities in a FRW cosmology 
is stable to small scalar, vector, and tensor inhomogeneities. We have shown 
that the existence of sudden singularities is stable when the density is bounded 
near the singularity except for some cases with special parameter choices. This 
result holds regardless of whether the background metric is spatially flat, closed 
or open. We also applied our analysis to a complementary characterization of 
Big Rip singularities and showed that they are stable if the leading term in the 
time-dependence of the FRW scale factor is proportional to a(t) = (t s — t) Vo 
where r] < —1, and unstable otherwise. As discussed in the introduction, 
there have been a number of identifications of sudden singularity occurrence 
in theories of gravity other than general relativity. The approach described in 
this paper can also be straightforwardly applied in these theories to determine 
whether the sudden singularities that occur there are also stable. 
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APPENDIX 
Equations of the form $ + AT S § + BT S <5> = 

Consider the equation 

$ + AT S <I> + BT S <5> = (64) 
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where $ = $(T), — 1 < s < and A, B are arbitrary constants. We want to 
investigate the behaviour of this equation as T — > 0. 

We first reduce the equation to canonical form <j> + f{T)4> = 0. The substi- 
tution <f> = </>cxp(— AT s+1 /2(s + 1)) achieves this, and yields 

BT S - ^AsT*- 1 - ^A 2 T 2s ^j = 0. (65) 
We now let Y(t) = t(p(l/t) and this gives 

Y + ^Bt~'- A - \Ast- s -* - ^A 2 t- 2s -^j Y = 0. (66) 

where t = 1/T and we now study the behaviour as t — > +oo. The substitutions 
Y(t) — exp(0(t)) and u(t) = <j>{t) transform (83) into a Riccati equation, which 
(after neglecting sub-leading terms) becomes: 

u + u 2 - ^Asr s - 3 = 0. (67) 

Finally we make the substitution u = y/y to obtain the Emden- Fowler equa- 
tion 

V = \Asyt- s -\ (68) 
This has a solution in terms of Bcsscl functions: 



(69) 



yi t) = «v» (c^ [z7?t' ( - s - 1)/2 J + [z^i' ( — 1 v 2 J J > 

(70) 

where D — As/2. We use the asymptotic behaviour of the Bessel functions as 
t — > +oo to obtain 

y(t) ~ C 3 + C 4 t (71) 
and it(t) - So Y(i) - t and $(T) -> consi as T -> 0. 
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